#Get objects
old_ws = ls()

#############
#Tables NVMS#
#############

useTable_nvms = mergedTable_nvms[!is.na(NVM_Violence_1_count) & (sample_ni_it)]

#################
#Events Analysis#
#################

#Model specifications
nb_specs = expand.grid(x = 'IT_win', y = c('NVM_Violence_1_count', 'NVM_Violence_1_deaths'), stringsAsFactors = F)
lm_specs = expand.grid(x = 'IT_win', y = c('NVM_Violence_1_binary', 'NVM_Violence_1_deaths_binary'), stringsAsFactors = F)


#Descriptives:
Map(function(y) summ_loop(y, useTable_nvms),
    c(nb_specs$y, lm_specs$y)) %>% 
  rbindlist %>%
  fwrite("./output/descriptives/nvms.csv")


#Estimate models
nvms_nb_models = Map(function(y,x) glmnb_loop(y,x, useTable_nvms),
                nb_specs$y, nb_specs$x)

nvms_lm_models = Map(function(y,x) lm_loop(y,x,  useTable_nvms),
                lm_specs$y, lm_specs$x)

nvms_list = c(nvms_nb_models, nvms_lm_models)
nvms_se = lapply(nvms_list, function(x) cluster_errors(x, useTable_nvms[, cluster]) %>% diag %>% sqrt)

########################
#Tables Violence: PODES#
########################

useTable_podes = mergedTable_podes[!is.na(Violence_1_count) & (sample_ni_it)]

#############
#Full sample#
#############

#Model specifications
nb_specs = expand.grid(x = 'IT_win', 
                       y = c(paste0('Violence_',1,'_count'),
                             paste0('Violence_',1,'_deaths')), stringsAsFactors = F
)

lm_specs = expand.grid(x = 'IT_win', 
                       y = c(paste0('Violence_',1,'_binary'), 
                             paste0('Violence_',1,'_deaths_binary')), stringsAsFactors = F
)


#Descriptives:
Map(function(y) summ_loop(y, useTable_podes),
    c(nb_specs$y, lm_specs$y)) %>% 
  rbindlist %>%
  fwrite("./output/descriptives/podes.csv")


#Estimate models
podes_nb_models = Map(function(y,x) glmnb_loop(y,x, useTable_podes),
                nb_specs$y, nb_specs$x)

podes_lm_models = Map(function(y,x) lm_loop(y,x,  useTable_podes),
                lm_specs$y, lm_specs$x)

podes_list = c(podes_nb_models, podes_lm_models)
podes_se = lapply(podes_list, function(x) cluster_errors(x, useTable_podes[, cluster]) %>% diag %>% sqrt)

############
#Make Table#
############

table_list = c(podes_list, nvms_list)
table_se = c(podes_se, nvms_se)

#create table
n_clusters = useTable_podes[(nvms_sample), cluster] %>% unique %>% length
n_clusters_2 = useTable_podes[, cluster] %>% unique %>% length

note_text = paste("Count models use negative binomial regression; binary outcomes use OLS. Standard errors are clustered by constituency-clusters. In the NVMS and PODES samples, there are ", n_clusters, "and ", n_clusters_2, " clusters, respectively. Observations are constituencies in which the last seat was contested by Islamist and secular nationalist parties with a margin less than 1 percent.")

table = stargazer(table_list, se = table_se, type = 'latex', 
                  title = "Estimated Effects of Islamist Victory on Religious Violence",
                  label = 'tab:main_violence_ni_it',
                  model.names = F,
                  model.numbers = T,
                  column.separate = c(4,4),
                  column.labels = c('PODES', 'NVMS'),
                  multicolumn = T,
                  dep.var.labels = rep(c("Events","Deaths"), 4), 
                  add.lines = list(c('Count', rep(c('Y', 'Y', 'N', 'N'), 2)),
                                   c('Binary', rep(c('N', 'N', 'Y', 'Y'), 2))),
                  covariate.labels = c("Islamist Win"),
                  star.cutoffs = c(0.1, 0.05, 0.01),
                  #float.env = 'sidewaystable',
                  notes.align = 'l',
                  font.size = 'scriptsize',
                  keep.stat = 'n',
                  style = 'apsr')



write_latex(table[-10] %>% append(table[10], after = 10) %>% append("\\cmidrule(lr){2-5}\\cmidrule(lr){6-9}", after = 10), note_text, './output/tables/table_1.tex')


#Clean up
drop = setdiff(ls(), c(old_ws)) 
rm(list = drop)